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Abstract: In this paper, we develop a random process theory to explain the 
laser speckle phenomena. The relation between the probability distribution 
of speckle's integrated intensity random process Y(f) and the relative 

velocity v t is derived. Based on the random process theory, traditional 

spatial or temporal laser speckle contrast analysis (i.e. spatial or temporal 
LASCA) can be derived as the spatial or temporal estimators respectively. 
Both spatial LASCA and temporal LASCA suffer from noise due to 
insufficient statistics and nonstationarity in either spatial or temporal 
domain. Furthermore, either LASCA results in a reduction of spatial or 
temporal resolution. A new random process estimator is proposed and able 
to overcome these drawbacks. In an in-vitro study, random process 
estimator outperforms either spatial LASCA or temporal LASCA by 
providing much higher SNR (random process estimator vs. spatial LASCA 
us. temporal LASCA: 33.64+6.87 ( mean±s.t.d. ) vs. 9.08+2.85 vs. 
3.83+1.05). In an in-vivo structural imaging study, random process 
estimator efficiently suppresses the noise in contrast image and thus 
improves the distinguishability of small vessels. In a functional imaging 
study of cerebral blood flow change in the somatosensory cortex induced by 
rat's hind paw stimulation, random process estimator provides much lower 
estimation errors in single trial data (random process estimator vs. temporal 
LASCA: 0.31+0.03 vs. 1.36+0.09) and finally leads to higher resolution 
spatiotemporal patterns of cerebral blood flow. 

©2009 Optical Society of America 

OCIS codes: (110.6150) Speckle imaging; (170.3880) Medical and biological imaging; 
(1 10.4280) Noise in imaging system. 
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1. Introduction 

Structural and functional brain imaging has become a very important tool in neuroscience 
research. As a simple two-dimensional imaging method, laser speckle imaging (LSI) [1] has 
been widely used in imaging the cerebral vasculature and cerebral blood flow (CBF) [2], e.g. 
the vasculature in brain tumor [3] and retina [4], the functional changes of CBF under focal 
cerebral ischemia [2], hypotension [5] and/or peripheral electrical stimulation [6]. Both 
spatial [7] and temporal [8] laser speckle contrast analysis methods, i.e. spatial LASCA and 
temporal LASCA respectively, have been used to estimate the contrast values ( K 2 ) which is 
inversely proportional to the relative velocity ( v) [9]. However, spatial LASCA decreases the 
spatial resolution while temporal LASCA is under the cost of temporal resolution. 
Furthermore, the inevitable noise in both LASCA methods not only decreases the 
distinguishability of vessels, but also results in errors in estimating functional changes of 
blood flow. 

In this paper, we develop a random process theory for explaining the speckle phenomena 
and derive the relation between the velocity of CBF and the probability distribution of 
speckle's random process. With the help of the random process theory, both LASCA methods 
are directly derived as the spatial and temporal estimators, and the noise is also modeled and 
analyzed. Finally, based on the property of noise, a new random process estimator is proposed 
to improve the signal-noise ratio (SNR) of the estimation, which thus leads to a better 
structural and functional imaging of CBF. 

2. Theory 

2.1 Background and preparation 

When coherent laser light illuminates a surface of rat's cortex, the random interference 
patterns, i.e. laser speckles, are produced by the coherent addition of scattered laser lights 
with slightly difference in light-path lengths. The motion of scattering particles, i.e. the red 
cells in blood vessels, changes the temporal pattern of speckles. The principle of LASCA is to 
estimate the relative velocity of CBF by analyzing the changing speckle patterns. 

Essentially, light is a kind of electromagnetic radiation. According to the electromagnetic 
theory, the area under coherent laser illumination can be modeled as a 'complex electric field' 
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E(t),t e[0,+<x>) [10,11]. Following the dynamic light scattering (DLS) theory [12], the 
scatters' velocity v(t) is related to E(t) 's autocorrelation function g^z) : 



g l (T) = (E(t)E\t-T)) l 



l{E{t)E\t)\ 



temporal ' 



(1) 



where r is the delay time; E*(t) is the complex conjugate of E(t) ; (•) temporal expresses the 
temporal average. The relation between v(f) and g, (r) relies on different models. For example, 
based on the assumption of Lorentzian distribution of velocity, Fercher and Briers [13] 
proposed the first model: g t (z) = exp(-r/r 0 ) (z 0 is the correlation time, r,, 1 <xv); based on 
diffusing wave spectroscopy (DWS) model [14], Bandyopadhyay et al. [15] suggested a more 
rigorous model: g,(r) = exp(— /(6t7 t q ) 112 ) (/is a constant parameter, z^ 1 ocv). In practice, 
it is very difficult to measure E(t) and g^z) directly, but it is common to study the 
intensity I(t),t e [0, +co) , defined as I{t) =\ E{t) I 2 . Generally, the recorded laser speckle 
image contains the intensity information of a composition of speckles with same speckle 
size s = 2.44/lf/# ( X is the wavelength of laser light, f/# is the / number of the imaging 

system) [2]. For each speckle, the scattering particles' velocity v(f) is related to£(f) and thus 
related to the statistical property of I(t) . 

2.2 Random process theory for laser speckle phenomenon 

Due to the randomness property of speckle phenomena, the intensity I(t), t e [0,+co) of one 
speckle is considered as a single realization of the continuous intensity random process X(f) . 
In this section, we are going to develop the random process theory to explain the laser speckle 
phenomenon and derive the relation between X(f) and velocity v(f) . 

At any time t, , I(t x ) is a single realization of the corresponding intensity random variable 
X(fj) whose probability distribution is related to the velocity v(f t ) . In practice, a CCD 
camera is used to acquire the intensity information of the imaging plane. By carefully 
adjusting the / number of the imaging system, we can make the pixel size of the recorded 

image equal to the speckle size so that there is only one speckle in one pixel. Therefore, the 
'intensity' value of any pixel is actually equal to the integrated intensity over the exposure 
duration T : 



where I(t) is a single realization of integrated intensity random process Y(f) = — I X(f ')df ' . 



Suppose the exposure time T is designed to be very short ( ~ 5ms ), then at any time f, , 
the velocity v(t'),t' e[t l ,t l +T] can be approximately considered as a constant. Thus, the 
distributions of random variables [X(t'),t' e [fpf, +T]} are the same, and each time-limited 
continuous random process X(f'), t ' e [t i , f, + T] can be considered as a wide-sense stationarity 
(WSS) random process. The properties of WSS random process [16] include: (a) its mean 
value ju x (f ') is a constant when t' e [f, , ^ + T] ; (b) its temporal average value of the product 

of any two samples, i.e. X(t') and X(t") , t', t" e [f,,f, +T] , depends only upon the time 
interval z = t"-t' . 

Based on the property (a), when the number of ensemble samples N — > +°o , we have: 




(2) 
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MO * <Y(f 1 )) mMe = l +T WW) mble = j { l+7 '(X( f ')) mWe dr' * Mx ( tl ) (3) 

Based on the property (b), the 'intensity autocorrelation function' g 2 (r) of X(f') during the 
exposure t' e [f,,f, +T] exists: 

g 2 (T) = (X(t')X(t' -r)> — / <X(f')>^ ora; . (4) 
According to the Siegert relation [12], we have: 

(X(t')X(t'-z)) temporal = (X(t l )) 2 enseM (l + /1\ gl (z) I 2 ), (5) 

where f' e [f t , f t + T] , 1/ /? is approximately the number of speckles per pixel (Here /? = 1 ). 
Suppose we obtain infinite temporal realizations of {X(t'),t' e[t l ,t l +T]} , then we have: 



<Y 2 (0> mWe =<^{ + ( + X{t')X{f)&t'dt") ememble 



(6) 

- (7 " 2 + f I> ^"-'') |2 df ' d O(/4(0>— 



<Y(0)L mHe (l + ^f2(l-f)^(r)dr) 



r * t 

Notes: (i) the temporal average (X(t')X(t")) lempoml in the time interval [t l ,t l +T] is used 
instead of the function X(t')X(t") in the second step based on the WSS property of 
{X(t'),t' e [f,,f, +T]} ; (ii) Eq. (5) is applied in the third step; (iii) Eq. (3) is used in the fifth 

step. 

By simply manipulating Eq. (6), we conclude: 



step; (iv) \ £ £/3\ gl (t"-t')\ 2 df'df" = ^|'2(l--)g 1 2 (r)dr [15,17] is used in the last 



ensemble 

<Y(^))L^ fe " ~T*> T gl 
Using Zun^^t,)).^ and Um N ^{Y\t l )) msM -(Y{t,))] nsM -^<7 2 (f,) , 

we rewrite Eq. (7) and define the contrast of integrated intensity random variable Y(f,) , i.e. 
K^{t x ) , as follows: 

Ki(t 1 ) = ^ = l[ ) 2(l-b g f(r)dr, (8) 

where g^r) is the £(?,) 's autocorrelation function [18]. 

Based on different models describing the relation between v(f) and g^f) , the relation 
between ^(f) and v(t) can also be derived directly [13,15,19,20]. For any single speckle, if 
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we can obtain its ensemble realizations (samples) {/ ( 2 (f,)}of the integrated intensity random 
variable Y(f t ) at time f, , the ensemble samples can be used to estimate the K\ and then the 
relative velocity v can be estimated. 

2.3 The spatial and temporal estimators and the estimation noise 

In practice, at any time?, , we can only get one realization of Y(f,) , i.e. I(t x ) , for each 
speckle (pixel). So we have to estimate K^it^ in other ways. An alternative method is to 
estimate the ju Y (f, ) and cr Y (f, ) using the intensities of the pixel ( /„ (t t ) ) and its surrounding 
pixels ( I j (t 1 ),j = l,...,M -1 ) based on the assumption that velocities at the surrounding 
pixels are very close. When this assumption is true, the distributions of integrated intensity 
variable {Y-(f,)} at different pixels are the same. Thus, spatial samples 

{Ij(t 1 ),j=0...M—l} can be considered as ensemble samples approximately: 
< Y (0)s Pa ,un=^Yt° I i *i *< Y ^>™^ and 
=m 2^o 7 , W*< Y • The estimations of // Y (f,) , <7 ¥ (0 and 

^y(^) can be obtained as follows: 

where / / s ( fl )=<Y(f 1 )> jpflfM *^ r (f 1 ), of (t x ) = <Y 2 (f, )) spatjal - <Y(f, )) 2 spatial «o$(f,). 

It is noted that this Eq. (9) is actually the spatial LASCA method [7]. Such estimation 
works well in most in-vitro simulation studies when the moving scattering particles in a large 
area have the same velocity. However, this assumption does not hold for a complex biological 
tissue, such as brain, where there is: 1) an abundance of scattering particles such as red cells, 
and 2) a large number of vessels with different diameters and velocities. 

Next, we analyze the estimation noise in spatial LASCA. Generally, the estimation of 
/u s {t{) and <t s (?j) in a spatial window centered at the interesting pixel could be affected by 

(i) limited size of the window, which results in the estimation noises NA^ (?,) and NA^ (f t ) ; 

and (ii) the spatial inhomogeneity of velocities within the window, which could result in 
estimation noises M^ (f,) and NB^it^ . Therefore, //,(?[) and cr^fj) can be modeled as 

follows: 

M^h) =MO+NA Ms (O +NB Ms (*,), (10) 

<r.<td = °V('i) + NA a s ^)+ NB a s (H) 
Although a larger window (e.g. 5x5 or larger) can result in smaller NA^ and NA t7s based 

on central limit theory, the NB^ and NB^ will increase due to higher inhomogeneity of 

velocities. In contrast, a smaller window (e.g. 3x3 ) will result in smaller NB„ and NB„ , but 

Ms u s 

larger NA^ and NA a ^ . 

Another problem of spatial LASCA is the loss of spatial resolution. When we need a high 
spatial resolution image of cerebral vessels, spatial LASCA is not appropriate. To preserve 
the spatial resolution, the temporal estimator was developed. 
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If we assume that the velocity in the interesting pixel is constant during the time 
interval ft,? 2 ] , then the probability distributions of integrated intensity variables 

[Y(t'),t' e ft,? 2 ]} are the same. In practice, we use a CCD camera to record the temporal 
samples ( {I(t')}) of integrated intensity random process Y(f') : 

1 rt,+nAt+T 

/ft+n A?) = -J[ +nA( I(t)dt', (12) 

where n = 0,...,N-1 , At is the time interval between two frames. Temporal samples 
{I(t f + nAt),n=0,...,N —1} are thus considered as ensemble samples of 

yc^wouw =^Z!: 0 l7 >. +« Af ) *< Y ( f .)>™* and 

<Y 2 ft)>,„, po „, =^Er/o 1/2 ^ +nAr ^< Y2 ^i)>™— • Then the estimation of iT Y ft) , i.e. 
^yft) and 0" Y ft) individually, is obtained by the temporal estimator defined as follows: 



Mi ft) 

where <u, ft ) = <Yft )>,^ ora/ « ^ ft ) , of ft ) = <Y 2 ft )) — - (Yft )>, 2 „ l?jora/ * a\ ft ) . 

Again the temporal estimator Kf empoml is actually the temporal LASCA method [8]. The 
temporal LASCA is also contaminated by the estimation noise, i.e. the noise components NA 
and NB in ju,{t { ) and <x,ft) as follows: 

A ft ) = lh ( f i ) + N \ (h ) + NB Ut ft ), (14) 

a-, ft ) = cr Y ft ) + N\ ft ) + Affl^ ft ). (15) 

where NA and A'A^, correspond to the noise due to limited statistics, NB and NB at 

correspond to the noise due to temporal inhomogeneity of velocities. 

Temporal LASCA preserves the spatial resolution because //,ft) and cr,ft) are 

calculated based on /ft + nAt)(n =0,...,N -I) of the same pixel. Therefore, high resolution 

structural imaging of cerebral blood vessel network is achieved. When imaging rat's CBF 
under stable condition, the velocity in each speckle is approximately stable during the 
recording, and the noise component NB can be ignored. However, if the blood velocity is 
changing, e.g. elicited in response to functional stimulation, the NB noise will be prominent. 

To achieve a high performance of statistical calculation (low NA noise), usually more 
than 50 frames are used in the temporal LASCA. In practice, the frame rate of the CCD 
camera is limited, for example IQfps in this study. So, the temporal resolution is 
compromised. Fewer frames, i.e. small time window, will reduce the loss of temporal 
resolution, but leads to larger NA noise. 

2.4 Random process estimator 

Based on the properties of integrated intensity random process Y(f) , we propose a new 

estimator, called random process estimator hereafter, which provides high SNR estimation 
and high spatiotemporal resolution for both structural and functional imaging. 
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The velocity of CBF can be considered as a constant within a short exposure duration (e.g. 
~ 5ms). For an interesting pixel, we calculate the {a s (t 1 +nAt),n =0,...,N -1} using a 3x3 
window from each frame. Because the window is small, there would be large noise 
component NA (Ts and small NB a in {cr s (f, +nAt),n =0,...,N-l} . Then, Eq. (11) is 

simplified into: <r s (fj + nAt) « <J Y (t l +nAt) +NA as (t l +nAt) , where n = 0,...,N -I. 

It is noted that each a Y (f, + nAt) is related to the velocity v(t l + nAt) and thus could 

change during the imaging procedure. The noise component NA a is due to the limited 

number of samples, so it is affected by the number of samples. Based on the central limit 
theorem, AM^ can be hypothesized to follow a zero-mean Gaussian distribution, i.e. 

N(0,af) , which will be confirmed based on the simulation data in Discussion. Therefore, the 
discrete sequences {NA^ {t l +nAt),n = 0,...,N -1} can be considered as independent and 

identically-distributed (IID) random variables with the same Gaussian distribution. 

Then we calculate the average of {a s (t l + nAt)} as the random process estimator for 
averaged <r Y (f, ) : 

Y N-l J JV-1 

<vd) = ttI>.A +wAf ) = °vd)+-I>Ax s ( f i +nt*), (16) 

" n=0 " „ =0 

where (Jyi^) = "^"o^y^i +nAt) corresponds to the averaged velocity v during 
[t l ,t l +(N -l)At] . Because [NA a (J, +nAt)} are IID random variables, based on the law of 

large number, lim — V^'AM (t, +nAt) — >0. So, there is a de-noising effect when N is 
greater than one. 

//,(f[) is used to estimate the averaged /7 Y ('i) m tne random process estimator, so there is 
no loss of spatial resolution in the calculation of ju (t^ . Actually, each /(f, +nAt) can be 
considered as a 'spatial estimation' of n s .(t l +nAt) with a lxl spatial window, which leads 
to a large noise component NA^ without NB^ noise. Ignoring the NB noise in Eq. (10), 

we have: /(;*, +nAt) = ju Y (t 1 +nAt) +NA JU ^ (t l +nAt) . Then we get: 

^ N-l ^ N-1 ^ N-1 

MveiO = T7Z^ f i +wAf ) = T7S^v(fi +"Ar) + — X^ a a (f i +nAt) - (17) 

A* „ = o Ar „ = o jv „ = o 

Similarly, we model NA^ i^ +nAt) as a sequence of IID random variables with a zero- 
mean Gaussian distribution based on the central limit theorem, i.e. NiO,^) (this assumption 
will also be confirmed based on the simulation data in Discussion). Therefore, when 

n^^o, ^Y.V a N \Mi +nM ^° and then ^( f i)=-"Y^)=^Z„ A '=<!/ / Y^+"Af) , 

where /7 Y ('i) corresponds to the averaged velocity v within [t lt t x + (N -l)Af] . 

Finally, the new random process estimator is used to estimate ^ Y (fj) for the averaged 
velocity v during time interval [fpfj + (N -l)At] as follows: 
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Compared with spatial LASCA, AT r 2 f provides much higher spatial resolution because the 
calculation of a uses only a 3x3 window. Furthermore, because the frame number N 

only affects the averaging (denoising) performance of the random process estimator, we can 
use fewer frames (e.g. A? =10) than the temporal LASCA (N = 50 images) to achieve a 
higher temporal resolution. Even using only 10 frames, the denoising performance of random 
process estimator is still significant while temporal LASCA is contaminated by both NA 
and NB . Therefore, random process estimator is very appropriate for functional imaging 
where both spatial and temporal resolutions need to be preserved as much as possible. 

Although temporal LASCA can provide high resolution image of blood vessel network 
based on a large number of raw images, the noise component NA is unavoidable. 
Furthermore, some physiologic changes, e.g. body temperature, heart rate, blood pressure and 
medicine, could also lead to a time-varying CBF and thus the NB noise, which will decrease 
the distinguishability of small vessels. On the other hand, larger N in random process 
estimator will have better denoising performance according to Eq. (16) and Eq. (17). 
Therefore, although random process estimator produces a slight loss in spatial resolution 
compared to temporal LASCA, it significantly improves the distinguishability of vessels. 

3. Experiments 

3.1 Imaging setup 

The imaging plane was illuminated with a 632nm He-Ne laser beam source (1508P-O, 
Uniphase, CA) which was reshaped by optical lens to expand the range of illumination. A 12- 
bit cooled CCD camera (Sensicam SVGA, Cooke, MI) with a 60mm macro (1:1 maximum 
reproduction ratio) / / 2.8 lens was held and focused on the imaging plane. Exposure time of 
the CCD camera was set to 5ms . Frame rate was 10 fps . The resolution of image recorded 
by the CCD camera was 704 x 704 pixels corresponding to an imaging area about 
4.7 x 4.7mm 2 . 

3.2 In-vitro simulation experiment 

A phantom is created to mimic blood vessels in vitro by flowing blood through plastic tubing 
using a syringe pump (PHD2000, Harvard Apparatus, MA). Real blood was obtained by ex- 
sanguination of rats destined for sacrifice and mixed with heparin sodium (10 IU/ml, Sigma 
Chemical Co., St. Louis, MO) to prevent clotting. A syringe filled with such blood was fitted 
into the syringe pump and connected to a polyethylene tube (internal diameter 0.034" , PE90). 
Blood was infused into the tube at 2mm/ s . When the flow was steady after ~ 5min , the tube 
was imaged using the imaging setup described above. 600 frames were acquired for analysis. 

3.3 In-vivo structural and functional LSI experiment 

The experimental protocol used in this study has been approved by the Animal Care and Use 
Committee of the Johns Hopkins Medical Institutions. The female Sprague-Dawley rats 
( ~ 320g ) were used in this study. To avoid the influence of surgical preparation on the rats' 
CBF, the thinned skull preparation was done on the first day. On the second day, the rats were 
taken back for the structural and functional LSI experiment. 

On the first day, the rat was anesthetized with intraperitoneal (IP) injection of mixture of 
90 mg/ kg of Ketamine and 10 mg/ kg of Xylazine. The animal was placed in a stereotactic 
frame (David Kopf Instruments, Tujunga, CA) throughout the experiment. The scalp was 
shaved and disinfected with 70% ethanol and povidone-iodine solution. All procedures were 
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performed using standard sterile precautions. After a midline scalp incision, the galea and 
periosteum overlying the parietal bone bilaterally were swept and retracted laterally. A 
5mmx5mm area centered at 2.5mm lateral and 2.5mm posterior to the bregma over the 
right cortex was thinned using a high speed dental drill (Fine Science Tools Inc. North 
Vancouver, Canada), until the inner cortical layer of bone was encountered. Rectal 
temperature was maintained at 37° C during the experiment using a homeothermic blanket 
system (model TP-500 T/Pump, Gaymar Industries, Inc., USA). 

On the second day, the rat was first held in a transparent chamber with 3% isofluorane 
gas and room air flow until becoming drowsiness. Its mouth and nose were then placed within 
a customized anesthesia mask in a well-fitting rodent size, which was connected to a C-Pram 
circuit designed to deliver and evacuate the gas through one tube. A mixed flow of 1.5% 
isofluorane, 80% oxygen and room air was delivered to the mask at the flow of 2L/min for 
anesthesia. 600 frames were recorded for the structural imaging experiment. 

Afterwards, subcutaneous needle electrode pairs (Safelead F-E3-48, Grass-Telefactor, 
West Warwick, RI) were inserted just under the skin of left hindpaw, one on the top of the 
hindpaw between 2 and 3 digits, the other one on the back skin of the hind paw, without 
direct contact to the nerve bundle. An isolated constant current stimulator (DS3, Digitimer 
Ltd., Hertfordshire, England) was used for the electrical stimulation of the hindpaw. Positive 
current pulses of 3.5mA magnitude and 200/xs duration at a frequency of 3Hz were used 

for hind paw stimulation. Each stimulation trial consisted of pre-stimulus 200 frames ( 20s), 
200 stimulation frames ( 20s ), and 400 post-stimulus frames ( 40s ), which was repeated ten 
times with several -minute in-between breaks. 

4. Results 

4.1 The denoising performance of random process estimator 

The denoising performance of random process estimator is evaluated in the in-vitro 
simulation experiment. During the experiment, the blood flow velocity in the polyethylene 
tube was strictly maintained at 2mm/s using the syringe pump. According to Eq. (16) and 

Eq. (17), the noise parts -^^J_ 0 ' NA U (t t + nAt) and ^ J_ 0 ' ^Ar s ( f i + n ^ ) m tne ran dom 

process estimator would be very closed to zeros with N = 600 . Therefore, for any pixel, 
H rpe and <r rpe based on all 600 raw images were considered as the true /u Y and cx Y . Figure 

1(a) shows one frame of the 600 raw images. Figure 1(c) is the estimated contrast image Ky 

using random process estimator. One pixel (the white point in the white box area in Fig. 1(c)) 
is selected arbitrarily to compare the estimation results based on spatial LASCA, temporal 
LASCA and random process estimator. 

For the selected pixel in Fig. 1(c), the result of spatial LASCA is calculated using a 7x7 
spatial window at each second (f, = 0,1,..., 595- ) while the results of temporal LASCA and 
random process estimator are calculated using a Is temporal window (10 frames) starting 
from each second. Figure 1(b) shows the estimations of cr Y based on spatial LASCA, 
temporal LASCA and random process estimator respectively. As expected, all estimations 
vary around the true value (white dashed line). Among all three methods, temporal LASCA 
produces the largest fluctuations because of large noise component NA (small temporal 
window) in Eq. (15). Although both LASCA methods are contaminated by NA noise, spatial 
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Fig. 1. The denoising performance of random process estimator: (a) one frame of raw images; 
(b) the estimation results of u v using different methods; (c) the contrast image K* estimated 

based on all 600 raw images using the random process estimator; (d) the estimated AT* using 
different methods. 



LASCA still provides a better estimation using a comparably bigger window size (spatial 
LASCA: TV = 7x7 = 49 ; temporal LASCA: A? = 10 ). Furthermore, according to the 
hydrodynamics, the blood velocity is still not distributed uniformly within the tube. 
Therefore, besides the NA noise, there is also NB noise in the result of both LASCA 
methods. Compared to the conventional LASCA methods, random process estimator provides 
the best estimation of a Y in the sense of smallest fluctuations. Since the estimations of /u Y 
are the same in temporal LASCA and random process estimator, we did not show the results 

of jU Y . 

Figure 1(d) shows the estimations of K Y at each second ( f, = 0,1,..., 59.? ) based on 

spatial LASCA, temporal LASCA and random process estimator where random process 
estimator outperformed the others again. To quantitatively compare the denoising 
performances of different methods, we define the signal-noise ratio (SNR) for the estimation 
of K Y as follows: 

SNR = (^) 2 = (— ^ ) 2 , (19) 

(n)-K 2 Y ff 2 

where T=60 in this study. For the selected pixel, K Y = 0.0053 , the SNRs are 2.75 (temporal 
LASCA), 10.34 (spatial LASCA), and 31.09 (random process estimator) respectively. For 
all 60x100 pixels in the white box area in Fig. 1(c), the averaged SNR of random process 
estimator is 33.64 + 6.87 (mean±s.t.d.) while the spatial LASCA and temporal LASCA are 
only with the averaged SNR 9.08 + 2.85 and 3.83 + 1.05 respectively. Clearly, random 
process estimator provides higher denoising performance than either spatial LASCA or 
temporal LASCA. 

Although spatial LASCA provides a slightly better estimation than temporal LASCA, the 
loss of spatial resolution is unavoidable and usually unacceptable. Since small vessels play an 
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important role in the in- vivo structural and functional studies, we'll only compare temporal 
LASCA and random process estimator in the rest sections. 

4.2 Random process estimator for structural imaging 

In the in-vivo experiment, 600 raw images ( 60s ) are recorded for imaging the structure of 
blood vessel network. During the imaging procedure, the rats were anesthetized and kept in 
stable condition, so blood velocities in the imaging area are assumed to be stable. Figure 2 
shows the contrast images estimated from the first 2 (Fig. 2(a,d)), 10 (Fig. 2(b,e)) and 80 (Fig. 
2(c,f)) frames using random process estimator and temporal LASCA respectively. For either 
method, the distinguishability of small vessels is improved when more frames are used. For 
example, the small vessel designated with black arrow in Fig. 2 (b) and (c) is hard to see in 
Fig. 2(a). 




Fig. 2. Structural imaging of rat's cerebral blood vessels: (a, b, c) show the contrast images 
estimated from the first 2, 10 and 80 frames using random process estimator; (d, e, f) 



show the contrast images K Y estimated from the first 2, 10 and 80 frames using temporal 
LASCA. 

According to Eq. (15), when more frames are used, the noise component NA decreases 
but still exists in the result of temporal LASCA. Furthermore, if there are small changes of 
blood velocity during the imaging procedure, the NB noise is also introduced in. However, 
using random process estimator, the noise part is rapidly averaged out according to Eq. (16), 
17) provided that more frames are used. More precisely, based on the same number of frames, 
random process estimator always provides better estimation and thus better distinguishability 
than temporal LASCA. For example, the small vessels in the white circled area in Fig. 2(c) 
are more distinguishable than in Fig. 2(f). 

The contrast values along the vertical white line in Fig. 2(c) estimated from different 
number of frames using the two methods is plotted in Fig. 3. There are six vessels denoted as 
V1-V6 crossing the line. In Fig. 3, the results of random process estimator are always 
smoother than temporal LASCA. In particular, based on only 2 frames, the noise in the result 
of temporal LASCA was prominent (Fig. 3(a)), while the result of random process estimator 
still reveals the denoising aspect especially in the vessel areas VI and V2. For 10 frames (Fig. 
3(b)), the random process estimator already provides an efficient estimation of K\ with high 
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Fig. 3. Denoising performance of random process estimator for structural imaging: the contrast 
values K 1 along the vertical line in Fig. 2 (c) estimated from different number of frames: N=2 
(a), N=10 (b), N=80 (c) using random process estimator and temporal LASCA. 
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Fig. 4. Spatiotemporal change of K during the pre-stimulation stage: the first row shows the 
K_ ,F ,...,F (from the left to the right) using random process estimator; the second raw 
shows the K , F F (from the left to the right) using temporal LASCA. 

SNR while the result of temporal LASCA is still very noisy. Furthermore, the denoising 
performance of random process estimator is better in the tissue areas compared to vessel 
areas. 

Although more frames can be used in the calculation by random process estimator, no 
significant improvement is found compared to Fig. 2(c). Therefore, 80 images are enough for 
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the structural imaging using random process estimator. In conclusion random process 
estimator provides better distinguishability of structural information of cerebral blood vessels. 
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Fig. 5. The averaged errors of single trial result using temporal LASCA and random process 
estimator respectively. 



4.3 Random process estimator for functional imaging 

In the functional stimulation experiment, each trial includes recording 800 raw images ( 80s ) 
continuously, i.e. 200 pre-stimulus frames ( 20s), 200 stimulation frames ( 20s) and 400 post- 
stimulus frames (40s). Images within each second (10 frames) are used to calculate one 
contrast image using random process estimator and temporal LASCA. Therefore, for each 
trial, we obtain 80 contrast images ( Kl 20 ,...,K^ 9 ). Kl 20 is used to normalize the change of 
contrast values as follows ( n = -19,... ,59): 



Because the stimulation procedure is 20s, other cortex areas may also be activated besides 
the somatosensory cortex. In this study, the area covering the activated somatosensory cortex 
is selected as the region of interest in the following discussions. 

To analyze the denoising performance for functional imaging, the {K^ 20 ,F_ 19 ,. . .,F_ ]b } of 

the interesting region in pre-stimulation stage of the first trial are shown in Fig. 4 using 
random process estimator (the first row) and temporal LASCA (the second row) respectively. 
Since the rat's condition was kept stable before the stimulation, there should not be 
significant changes in the contrast images and thus the values in images [F n ,n = -19, ...,-1} 

are expected to be much closed to zero. However, as shown in Fig. 4, temporal LASCA leads 
more noises than random process estimator. To quantify this estimation error, we define the 
'averaged errors' in each F n as follows: 




i X 100%, 



(20) 



1 



<V„ "col 



Error F =( 



(21) 



row" col x=l y=l 
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Fig. 7. The changing patterns of CBF under functional stimulation calculated by temporal 
LASCA. 

Figure 5 shows the averaged errors of single trial estimation using temporal LASCA and 
random process estimator respectively. For random process estimator, the averaged errors in 
all {F,n = -19,...,-1} is only 0.31 + 0.03 (mean±sX.d.), in contrast with 1.36±0.09 for 
temporal LASCA. 



#1 19045 - $15.00 USD Received 26 Oct 2009; revised 7 Dec 2009; accepted 14 Dec 2009; published 23 Dec 2009 
(C) 2010 OSA 4 January 2010 / Vol. 18, No. 1 / OPTICS EXPRESS 232 



{F Q ,...,F 59 } correspond to functional changes during the stimulation and post-stimulation 
procedure. By averaging the {F Q ,...,F 59 } of all ten trials, the averaged spatiotemporal 
changing map ({F 0 ,...,F S9 } ) is obtained. Figure 6 and Fig. 7 show the averaged changing 
map{F 0 ,...,F 33 } , using random process estimator and temporal LASCA respectively. As a 

reference, the first frames in Fig. 6 and Fig. 7 are the structural image calculated by random 
process estimator. 

In the functional imaging experiment, the blood flow in the somatosensory cortex changes 
in response to electrical stimulation. Therefore, besides the NA noise, the NB noise is also 
introduced into the result of temporal LASCA. Hence, the result of temporal LASCA (Fig. 7) 
is too noisy to precisely quantify the changes of CBF in small vessels. However, random 
process estimator (Fig. 6) provides high SNR and high spatiotemporal resolution for the 
functional changing pattern of CBF. 
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Fig. 8. The functional changes of contrast values corresponding to the pixel designated by a 
white circle in the first frame of Fig. 6: (a) the averaged changes, i.e. F ,n = -19, ...,59s , 
using temporal LASCA and random process estimator; (b) the standard deviation of the 
F ,n = -19, ...,59 s across all ten trials using temporal LASCA and random process 
estimator. 

In Fig. 6, the response of CBF induced by functional stimulation starts from 4s and 
reaches maximum at 20s ~ 22s , then recovers to the baseline at 32s . All significant changes 
can be located to small vessels (e.g. V1-V4 in Fig. 6). VI is the vessel with the most 
significant changes, which is activated earlier than other vessels and shows longer recovery 
duration. With the help of high spatial resolution, we also find that not all vessels are 
activated in this area (e.g. V5). 

One pixel in vessel VI is chosen for quantitative analysis (the white circle in the first 
frame of Fig. 6). The averaged value and standard deviation of {F n ,n = -19,. ..,59s} at this 
pixel across all 10 trials are plotted in Fig. 8 for both methods. The result of temporal LASCA 

exhibits a noisy fluctuations ( -10% 40% ) even in the pre-stimulation stage, so that only 

changes greater than 50% can be considered as significant change during stimulation. 
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However, the result of random process estimator shows much smaller fluctuations and a more 
stable profile. Furthermore, the standard deviation of the {F n ,n = -19,... ,59s} over all 10 
trials using random process estimator are always less than that using temporal LASCA (Fig. 8 
(b)), implying that random process estimator also suppress the noisy variations over trials. 




Fig. 9. Gaussian assumptions confirmed by the data acquired in the in-vitro experiment: (a) the 
pdf of NA„ ; (b) the P-P plot of all 600 NA„ s; (c) the pdf of NA n ; (d) the P-P plot of all 

600 NA„ s. 
5. Discussion 

5.7 The Gaussian noise assumptions for random process estimator 

In the development of random process estimator, the most important hypotheses are the 
Gaussian noise models in Eq. (16) and Eq. (17). Using the data acquired in the in-vitro 
experiment, we could test the hypotheses. 

For the selected pixel in Fig. 1(c), 600 estimations of tr Y can be obtained using a 3x3 

spatial window. The noise component NA Us of each estimation is then calculated by 

subtracting the true value cx Y . The probability density function (pdf) is then calculated from 

all 600 NA t7s by the kernel smoothing density estimation (Fig. 9(a)). The pdf is very close to 

a Gaussian distribution with 0±2.158xl0" 3 ( mean + s.t.d. ). The P-P plot of all 600 NA^ in 

Fig. 9(b) is also very close to that of a Gaussian distribution. The hypothesis that the 
distribution of NA Us is a Gaussian distribution is further confirmed by the K-S test with a 

high z value ( z = 0.859 ). 

Similarly, the distribution of the noise component NA^ , i.e. I. — fj^ (i = 1, . . ., 600) , is also 

tested by K-S test with a high z = 0.848 . The pdf of the distribution ( 0 + 9.405 x 10" 3 ) and the 
P-P plot are shown in Fig. 9(c) and (d) respectively. These tests confirm that the distribution 
of NA,. is also a Gaussian distribution. 

Ms 
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5.2 Comparison with hybrid temporal and spatial method 

Recently, a hybrid temporal and spatial method was proposed to improve the SNR by 
spatially averaging (5x5 spatial window) the result of temporal LASCA which was 
successfully applied in imaging rat's retinal blood flow [21]. Similar to spatial and temporal 
LASCA, the hybrid method works well when blood flow is constant within both spatial and 
temporal window, i.e. NB noise does not exist. We further process the same data set in our 
study using the algorithm in [21] to compare the performances of random process estimator 
and hybrid method in SNR, distinguishability of small vessels and revealing the functional 
responses in small vessels. 

In the in-vitro simulation experiment, for the selected pixel in Fig. 1(c), the hybrid method 
produces a higher SNR (16.48) than temporal LASCA (2.75) and spatial LASCA (10.34), 
which, however, is still significantly lower than the SNR by random process estimator 
(31.09), because the use of a 5x5 spatial window increases the NB noise while suppressing 
the NA noise. In the structural imaging of the selected area in Fig. 10(a), the details of small 
vessels (indicated with black arrows) would be hardly visible by the hybrid method with 
a5x5 spatial window (Fig. 10(b)). Although using a 3x3 spatial window in hybrid method 
slightly improves the spatial resolution (Fig. 10(c)), but the distinguishability of the vessels is 
still not as good as that by random process estimator (Fig. 10(d)) due to the decrease of SNR. 
Furthermore, the hybrid method cannot suppress the NB noise when there are significant 
blood flow changes within the temporal window. For example, NB noise is so large that the 
small functional changes are strongly contaminated by background noise in temporal LASCA 
(Fig. 1 1(b)) under the functional stimulation. Hybrid method will filter out noise and the weak 
functional signals as well (Fig. 11(c)). However, random process estimator well reveals the 
functional changes in the small vessel while removing the noise (Fig. 1 1(d)). 




Fig. 10. The vascular details in the white box area in (a) using hybrid method with 5x5 (b) 
and 3x3 (c) spatial window, random process estimator (d). 




Fig. 11. The functional CBF changes of selected area in (a) at the 8th sec estimated by 
temporal LASCA (b) as in Fig. 7, hybrid method (c) and random process estimator (d) as in 
Fig. 6. 



6. Conclusion 

In this paper, we develop the random process theory for the laser speckle phenomenon and 
thus propose a new random process estimator to estimate the true value of . Compared 
with traditional spatial LASCA and temporal LASCA, random process estimator provides 
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higher SNR and higher spatiotemporal resolution contrast images for both structural and 
functional imaging of cerebral blood vessels and flow. 
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